
#########################
## Load and clean data ##
#########################

# load data
trdmat <- read.csv("data/datasets/main_data.csv")
excl_agg <- read.csv("data/datasets/excl_agg.csv")
tars_covered <- read.csv("data/datasets/tars_covered.csv")

# load export-based measure 
exp <- read.csv("data/datasets/expoth.csv")
exp$hs6code <- as.character(exp$hs6code)
exp$hs6code[nchar(exp$hs6code) == 5] <- paste("0", exp$hs6code[nchar(exp$hs6code) == 5], sep = "")

# create a version of trdmat that is aggregated up to the right level
trdmat_hs6 <- aggregate(cbind(impchn20122016,impchn20182020,impchn20212023,impALL20122016,impALL20182020,impALL20212023,impwor20122016,impwor20182020,impwor20212023,
  impALL20122016allies,impALL20182020allies,impALL20212023allies,impALL20122016nonallies,impALL20182020nonallies,impALL20212023nonallies,
  impALL20122016friends,impALL20182020friends,impALL20212023friends,impALL20122016nonfriends,impALL20182020nonfriends,impALL20212023nonfriends) ~ 
  substr(hs10code, 1, 6), FUN = sum, data = trdmat)
names(trdmat_hs6)[1] <- "hs6code"
trdmat_hs6 <- join(trdmat_hs6, exp, type = "left")

######################################
## Models imports 20122016-20182020 ##
######################################

## top 10
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20182020/3, trdmat_hs6$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*market*I(impothpos/10), data = newdat2);
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top10 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20182020/3, trdmat_hs6$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2); 
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top25 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20182020"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20182020/3, trdmat_hs6$impALL20182020/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2);
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALL <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs6code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top10[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_top10_20182020 <- tab

hs6codes <- substr(newdat2_top25$hs6code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top25[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_top25_20182020 <- tab

hs6codes <- substr(newdat2_ALL$hs6code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_ALL_20182020 <- tab

######################################
## Models imports 20122016-20212023 ##
######################################

## top 10
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20212023/3, trdmat_hs6$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016,trdmat_hs6$exptop10pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_top10 <- newdat2
mod <- lm(log(imp+1) ~ year*market*I(impothpos/10), data = newdat2);
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top10 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2top10 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## top 25
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20212023/3, trdmat_hs6$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016, trdmat_hs6$exptop25pos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_top25 <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2);
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1top25 <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2top25 <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## ALL
newdat <- data.frame(year = rep(c("20122016","20212023"), each = 2*nrow(trdmat_hs6)), 
        hs6code = trdmat_hs6$hs6code,
        imp = c(trdmat_hs6$impchn20122016/5, trdmat_hs6$impALL20122016/5, trdmat_hs6$impchn20212023/3, trdmat_hs6$impALL20212023/3)/1000,
        market = c(rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6)), rep("China", nrow(trdmat_hs6)), rep("Other", nrow(trdmat_hs6))),
        impothpos = c(trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016,trdmat_hs6$expALLpos20122016)
)
newdat$market <- factor(newdat$market, levels = c("Other","China"))

postrade <- trdmat_hs6$hs6code[trdmat_hs6$impchn20122016/5 > 0 | trdmat_hs6$impwor20122016/5 > 0]
newdat2 <- na.omit(newdat[newdat$hs6code %in% postrade,])
newdat2_ALL <- newdat2
mod <- lm(log(imp+1) ~ year*market*log(impothpos+1), data = newdat2);
mm_df <- data.frame(model.matrix(mod), imp = log(newdat2$imp+1), hs6code = newdat2$hs6code)
mm_df <- na.omit(mm_df)
colnames(mm_df) <- c("Intercept","year","market","impothpos",
  "year_market","year_impothpos","market_impothpos","year_market_impothpos",
  "imp","hs6code")
aggvars <- apply(mm_df[,colnames(mm_df)[2:9]], 2, as.numeric)
mm_df$hs6code_market <- paste(mm_df$hs6code,mm_df$market, sep = "")
agg <- aggregate(aggvars ~ hs6code_market, FUN = mean, data = mm_df)
colnames(agg) <- c("hs6code_market","year_mean","market_mean","impothpos_mean",
  "year_market_mean","year_impothpos_mean","market_impothpos_mean","year_market_impothpos_mean",
  "imp_mean")
mm_df <- join(mm_df, agg)
mm_df[,2:9] <- mm_df[,2:9] - mm_df[,12:19] 

mod1ALL <- lm(imp ~ year + year_market + year_impothpos + year_market_impothpos, data = mm_df)
mm_df$year_sector <- paste(newdat2$market, newdat2$year, substr(newdat2$hs6code,1,2),sep = "_")
mod2ALL <- lm(imp ~ year_impothpos + year_market_impothpos + year_sector, data= mm_df); 

## write models
hs6codes <- substr(newdat2_top10$hs6code, 1, 6)
mods_top10 <- list(coeftest(mod1top10, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top10, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top10[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top10), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top10)){mod <- mods_top10[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top10)$df[1:2]),sum(summary(mod2top10)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_top10_20212023 <- tab

hs6codes <- substr(newdat2_top25$hs6code, 1, 6)
mods_top25 <- list(coeftest(mod1top25, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2top25, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_top25[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_top25), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_top25)){mod <- mods_top25[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1top25)$df[1:2]),sum(summary(mod2top25)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_top25_20212023 <- tab

hs6codes <- substr(newdat2_ALL$hs6code, 1, 6)
mods_ALL <- list(coeftest(mod1ALL, vcov = vcovCL, cluster = hs6codes)[1:5,1:4], 
  coeftest(mod2ALL, vcov = vcovCL, cluster = hs6codes)[1:7,1:4])

vars <- c(rownames(mods_ALL[[1]])[c(5)])
tab <- matrix(data = NA, ncol = length(mods_ALL), nrow = length(vars)*2)
rownames(tab) <- c(rep(vars, each = 2))
for(i in 1:length(mods_ALL)){mod <- mods_ALL[[i]]; tabinmod <- unique(rownames(tab)[rownames(tab) %in% rownames(mod)])
stars <- rep("", length(tabinmod)); stars[mod[tabinmod,4] < .1] <- "^{+}"; stars[mod[tabinmod,4] < .05] <- "^{*}"; 
stars[mod[tabinmod,4] < .01] <- "^{**}"; stars[mod[tabinmod,4] < .001] <- "^{***}";  
tab[rownames(tab) %in% rownames(mod),i] <- c(rbind(paste(myround(mod[tabinmod,1],3), stars, sep = ""), paste("(", myround(mod[tabinmod,2],3), ")", sep = "")))}
rownames(tab) <- c(rbind(vars, "")); rownames(tab)[c(1)] <- 
  c("Post-Trade War$\\cdot$China$\\cdot$Alternative markets") 
N <- c(sum(summary(mod1ALL)$df[1:2]),sum(summary(mod2ALL)$df[1:2]))
N <- paste("\\multicolumn{1}{c}{", N, "}", sep = ""); tab <- rbind(tab, N)
addtorow <- list(); addtorow$pos <- list(2,3); addtorow$command <- c(paste0('\\midrule '),paste0(''))

tab1_ALL_20212023 <- tab

ptable(cbind(rownames(tab), tab1_top10_20182020, tab1_top10_20212023), "paper/tables/TableA2top.tex")
ptable(cbind(rownames(tab), tab1_top25_20182020, tab1_top25_20212023), "paper/tables/TableA2mid.tex")
ptable(cbind(rownames(tab), tab1_ALL_20182020, tab1_ALL_20212023), "paper/tables/TableA2bot.tex")

############
## Figure ## 
############

# Figure 2
pdf("paper/figs/FigureA1.pdf", height = 6, width = 3.5)
par(mfrow = c(3,1), mar=c(0,1.5,1.25,.5), oma = c(.5, .2, 0, 0), lwd = .5) # bottom, left, top, right

barplot(table(trdmat_hs6$exptop10pos20122016[trdmat_hs6$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-100,1500))
text(x = 1.2*seq(0,9,1) + .7, y = -35, labels = seq(0,9,1), cex = .6)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6, labels = seq(0,1500,500), at = seq(0,1500,500))
title(main="# of non-China top-10 markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

barplot(table(trdmat_hs6$exptop25pos20122016[trdmat$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-10,150))
text(x = 1.2*seq(0,24,1) + .7, y = -5, labels = seq(0,24,1), cex = .5)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6)
title(main="# of non-China top-25 markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

barplot(table(trdmat$impALLpos20122016[trdmat$impchn20122016 > 0]), axes = F, main = "",
  names.arg = "", ylim = c(-80,685), xlim = c(0,200), border = "gray")

text(x = 1.2*seq(0,160,10) + .7, y = -47, labels = seq(0,160,10), cex = .5)
axis(2, cex.axis = .4, tck = -.015, padj = 3.6, labels = seq(0,600,100), at = seq(0,600,100))
title(main="# of non-China markets producing product", cex.main = .6, line = -.5)
mtext("Number of products", 2, line=.85, cex = .5)

dev.off()

